function [evals, evecs] = slsymeig(A, k, ord, varargin)
%SLSYMEIG Compute the eigenvalues and eigenvectors for symmetric matrix
%
% [ Syntax ]
%   - evals = slsymeig(A)
%   - [evals, evecs] = slsymeig(A)
%   - ... = slsymeig(A, k)
%   - ... = slsymeig(A, k, ord, ...)
%
% [ Arguments ]
%   - A:        the target symmetrix matrix to be decomposed
%   - k:        the number of eigenvalues to be solved.
%               If k is empty or omitted, it is set to the dimension
%               of the matrix.
%   - ord:      the order of the eigenvalues in the result 
%               \{:
%                   - descend: the eigenvalues are sorted from large
%                              to small
%                   - ascend:  the eigenvalues are sorted from small
%                              to large
%               \:}
%               If ord is empty or omitted, by default, it is set to
%               'descend'.
%   - evals:    the column vector of eigenvalues of A.               
%   - evecs:    the matrix of eigenvectors of A
%               the eigenvectors are arranged in the order corresponding
%               to that of eigenvalues.
%
% [ Description ]
%   - evals = slsymeig(A) computes only the eigenvalues of symmetric 
%     matrix A.
%
%   - [evals, evecs] = slsymeig(A) computes both the eigenvalues and
%     eigenvectors of matrix A. If A is a d x d matrix, then evals
%     is a d x 1 column vector with the eigenvalues in descending order.
%     evecs is a d x d matrix with each column vector corresponding to
%     an eigenvalue in evals.
%
%   - ... = slsymeig(A, k) the user can specify the number of eigenvalues
%     to preserve, if not specified, all eigenvalues will be preserved.
%
%   - ... = slsymeig(A, k, ord, ...) the user can specify how to sort the
%     eigenvalues. By default, they are sorted in descending order.
%
%     You can further specify more options to control the solution
%     \{:
%         - use_eigs:     whether to use eigs to solve eigen-problem
%                         By default, it is []. The function chooses to
%                         use_eigs when k <= d/5 or A is sparse.
%                         You may also set it to true or false.
%         - force_sym:    whether to enforce that A is symmetric before
%                         solving. (default = true)
%                         The program enforces this by 0.5 * (A + A')         
%     \:}
%
% [ History ]
%   - Created by Dahua Lin on Apr 21, 2006
%   - Modified by Dahua Lin on Sep 9, 2006
%       - The user now can specify the number of eigenvalues to preserve
%         if necessary.  
%       - An auto-selection scheme, which will use eigs when k is small
%         in order to achieve higher efficiency.
%   - Modified by Dahua Lin on Jul 15, 2007
%       - fix the help
%       - change the argument checking   
%       - refactor the program structure.
%       - add options: use_eigs and force_sym.
%

%% parse and verify input arguments

assert(isfloat(A) && ndims(A) == 2, 'sltoolbox:slsymeig:invalidarg', ...
    'A should be a symmetric matrix.');
[d, d2] = size(A);
assert(d == d2, 'sltoolbox:slsymeig:invalidarg', ...
    'A should be a symmetric matrix.');

if nargin < 2 || isempty(k)
    k = d;
else
    assert(isscalar(k) && k > 0 && fix(k) == k, 'sltoolbox:slsymeig:invalidarg', ...
        'k should be empty or a positive integer scalar.');
    assert(k <= d, 'sltoolbox:slsymeig', ...
        'k exceeds its upper bound %d', d);
end


if nargin < 3 || isempty(ord)
    ord = 'descend';
else
    assert(any(strcmp(ord, {'descend', 'ascend'})), 'sltoolbox:slsymeig:invalidarg', ...
        'ord should be either ''descend'' or ''ascend''.');
end

opts = struct( ...
    'use_eigs', [], ...
    'force_sym', true);
opts = setopts(opts, varargin{:});

assert(isempty(opts.use_eigs) || (islogical(opts.use_eigs) && isscalar(opts.use_eigs)), ...
    'sltoolbox:slsymeig:invalidopt', 'The option use_eigs should be empty or a logical scalar.');
if isempty(opts.use_eigs)
    use_eigs = (k <= d / 5 || issparse(A));
else
    use_eigs = opts.use_eigs;
end

assert(islogical(opts.force_sym) && isscalar(opts.force_sym), ...
    'sltoolbox:slsymeig:invalidopt', 'The option force_sym should be a logical scalar.');

%% compute

% enforce symmetry
if opts.force_sym
    A = 0.5 * (A + A');
end

% decide scheme and compute
if ~use_eigs
    [evecs, evals] = eig(A);
else
    [sigma, eopts] = get_eigs_opts(ord);
    [evecs, evals] = eigs(A, k, sigma, eopts);
end
    
evals = diag(evals);

% enforce real
if ~isreal(evecs)
    evecs = real(evecs);
end
if ~isreal(evals)
    evals = real(evals);
end


%% sort: make the eigenvalues in descending order

[evals, si] = sort(evals, 1, ord);
evecs = evecs(:, si);

k0 = length(evals);
if k < k0
    evals = evals(1:k);
    evecs = evecs(:, 1:k);
end


%% Auxiliary function

function [sigma, eopts] = get_eigs_opts(ord)

switch ord
    case 'descend'
        sigma = 'LM';
    case 'ascend'
        sigma = 'SM';
end

eopts = struct('disp', 0, 'issym', true);





